function obj = calibrate_lapse_u(lapse_u_input, D_input, N,...
    alpha, ...
    Delta_xi_all,...
    K, pr_y, N_y, resource_y1, resource_y2, y_medi, x_fc, pr_fc, ...
    rho, crra, sigma_c1, beta_c,  T1, B1_c, B2_c, delta, ...
    myopic, id_ist_old)
%% use data for prices and fringe firms
p1_ini = D_input.p1; % NN x 1
p2_ini = D_input.p2; % NN x KM
n_ini_unmerge = D_input.n_insurers; % NN x 1: 1 for majors, # fringes for the fringe
n_ini = zeros(N.st,1); % N_st x 1
for mm=1:N.st
    n_ini(mm) = unique(n_ini_unmerge(D_input.mkt==mm & D_input.major==0));
end

%% lapse+edit: codes below copied from solve_eq.m and then edited

%% notations
% N.ist: # majors + # markets 
% N_ist: # majors + # fringes (unmerged version)

%% consumer's value for outside option
state_pr_s = unique([D_input.state D_input.pr_k], 'rows'); 
% lapse+edit
lapse_u = lapse_u_input; % set to the current value 

for ss = 1:N.s

    pr_s = state_pr_s(ss,2:K+1)';  % pr_s: K x 1

    % consumer_v0_ftn called from demand estimation folder
    [v0, U0] = consumer_v0_ftn(N_y, rho, crra, beta_c, K, ...
        pr_s, x_fc, pr_fc, y_medi, B1_c, B2_c, T1, alpha,...
        lapse_u, resource_y1, resource_y2);

    S(ss).v0 = v0; % N_y x 1
    S(ss).U0 = U0; % N_y x K
end

%% solve the model for each market
KM = size(D_input.p2,2);

max_p = min([resource_y1; resource_y2])-1;

% loop through each market (calendar year x geographical state combination)
lapse_model_mkt = zeros(N.st,1);

% store just to check 
s1_mkt_all = zeros(N.ist, 1);

for mm = 1:N.st
    %% indices for firms in market mm
    % index for majors
    iind_m = find(D_input.mkt==mm & D_input.major==1); % iind_m, = 1,..., N_ist
    NN_m = length(iind_m);

    % index for fringes
    iind_f = find(D_input.mkt==mm & D_input.major==0); % iind_f = 1,..., N_ist

    % needed to store in N.ist arrays with the original sorting as D
    id_ist_m = D_input.id_ist(iind_m); % NN_m x 1
    id_ist_f = unique(D_input.id_ist(iind_f)); % 1 x 1
    if length(id_ist_f)~=1
        disp('EROR in length(id_ist_f)')
    end

    %% exogenous values
    state_mm  = unique(D_input.state(iind_m)); % geographical state

    pr_s_row = D_input.pr_s(iind_m(1),:); % 1 x KM: the prob dist of agg states

    Delta_xi_m = Delta_xi_all(iind_m,:);
    Delta_xi_f = Delta_xi_all(iind_f(1),:);

    % outside option
    v0_mm = S(state_mm).v0;  % N_y x 1
    U0_mm = S(state_mm).U0;  % N_y x K

    %% initial guess on optimal choices
    p1_m = p1_ini(iind_m,:); % NN_m x 1
    p1_f = p1_ini(iind_f(1),:); % just one row for the symmetric fringe
    p2_m = p2_ini(iind_m,:); % NN_m x KM
    p2_f = p2_ini(iind_f(1),:); % just one row for the symmetric fringe

    n_variety = n_ini(mm); % 1 x 1, # fringes in the market

    p1_m(p1_m>max_p)=max_p;
    p1_f(p1_f>max_p)=max_p;
    p2_m(p2_m>max_p)=max_p;
    p2_f(p2_f>max_p)=max_p;

    %% demand using (p1, p2, n_variety) from the previous iteration
    NN = NN_m + n_variety;

    % stack majors followed by repeated fringes
    p1 = [p1_m; repmat(p1_f,n_variety,1)]; % NN x 1
    p2 = [p2_m; repmat(p2_f,n_variety,1)]; % NN x KM
    Delta_xi = [Delta_xi_m; repmat(Delta_xi_f, n_variety,1)];

    % correct beliefs
    % - v_c: N_y x NN
    p2_input = p2; % NN x KM
    [v_c, lapse_NN_c, ~] = predict_v_lapse(p1, p2_input, ...
        alpha,  Delta_xi,...
        N_y, resource_y1, resource_y2, pr_s_row, KM,...
        rho, crra, beta_c, U0_mm, B1_c, B2_c, T1, delta);

    % misinformed beliefs
    % - v_m: N_y x NN
    p2_input = repmat(p1, 1, KM); % NN x KM
    [v_m, ~, ~] = predict_v_lapse(p1, p2_input, ...
        alpha,  Delta_xi,...
        N_y, resource_y1, resource_y2, pr_s_row, KM,...
        rho, crra, beta_c, U0_mm, B1_c, B2_c, T1, delta);

    % choice probabilities including outside option -> per-fringe
    % shares are computed
    ind_temp = ones(1,NN);
    s1_ind_c = share_ind_ftn(v_c, ind_temp, v0_mm, sigma_c1); %N_y x NN+1
    s1_ind_m = share_ind_ftn(v_m, ind_temp, v0_mm, sigma_c1); %N_y x NN+1

    % N_y x NN+1
    s1_ind = (1-myopic)*s1_ind_c + myopic*s1_ind_m;

    s1_ind1 = s1_ind(:,1:end-1); % N_y x NN (inside options only)
    s1_mkt = sum(s1_ind1.*repmat(pr_y, 1, NN))'; % NN x 1, integrated over y

    % predict lapse probability
    % even consumers who were misinformed ex-ante will make lapse decisions
    % based on actual price increase. 
    lapse_NN = lapse_NN_c; 
    % - type distribution among the insured for each firm
    weight_y = s1_ind1./sum(s1_ind1, 1); % N_y x NN
    % - for each firm, average lapse prob using type distribution among the
    % insured
    lapse1 = sum(lapse_NN.*weight_y,1)'; % NN x 1
    % - average firm-level lapse prob using mkt share as the weight
    weight_firm = s1_mkt./sum(s1_mkt); % NN x 1
    lapse_model_mkt(mm) = sum(lapse1.*weight_firm);

    % for debugging 
    s1_mkt_m = s1_mkt(1:NN_m,:);
    s1_mkt_f = s1_mkt(NN_m+1,:);
    for jj=1:NN_m
        ll=find(id_ist_m(jj)==id_ist_old);
        s1_mkt_all(ll,:) = s1_mkt_m(jj,:);
    end
    ll = find(id_ist_f==id_ist_old);
    s1_mkt_all(ll,:) = s1_mkt_f;
end

%% evaluate the difference: predicted vs. empirical lapse rate
obj = (mean(lapse_model_mkt) - unique(delta))^2;


















